This dataset contains many different animals across space and time, recording a variety of data that seems meaningful to the groups. In this project I chose to work with primates.

The primates in this dataset contain Adult Body Mass, and Adult Head-to-Body lengths, which is what I have chose to model, along with using their taxonomic groups. The hopes are to model the relationship of the animal’s size(lengths) to their mass, while considering their groups in a multilevel model.

Load Libraries

library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.12.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:stats':
## 
##     ar
library(dagitty)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(reshape2)
library(ggplot2)
library(RColorBrewer)

#Load Dataset

df.animals <- read.csv2("data/BoneLengths.csv", sep=",", fileEncoding="UTF-8-BOM",as.is=T) #Encoding fixes name ERR
df.animals
#What order has the most data?
#counts <- df.primate %>% group_by(df.primate$MSW05_Order) %>% count() %>% arrange(desc(n))

counts <- df.animals %>% group_by(df.animals$MSW05_Genus) %>% count() %>% arrange(desc(n))
counts
df.animals$X5.1_AdultBodyMass_g <- as.numeric(df.animals$X5.1_AdultBodyMass_g)
df.animals$X8.1_AdultForearmLen_mm <- as.numeric(df.animals$X8.1_AdultForearmLen_mm)
df.animals$X13.1_AdultHeadBodyLen_mm <- as.numeric(df.animals$X13.1_AdultHeadBodyLen_mm)

Primates were chosen from the dataset

df.animals <- df.animals %>% filter(df.animals$X5.1_AdultBodyMass_g != -999)
df.animals <- df.animals %>% filter (df.animals$X5.1_AdultBodyMass_g != 1)
#df.animals <- df.animals %>% filter (df.animals$X8.1_AdultForearmLen_mm != 1)
df.animals <- df.animals %>% filter (df.animals$X13.1_AdultHeadBodyLen_mm != 1)
df.animals <- df.animals %>% filter (df.animals$X13.1_AdultHeadBodyLen_mm != -999)

df.primate <- df.animals %>% filter(df.animals$MSW05_Order == "Primates")
#df.bats <- df.animals %>% filter(df.animals$MSW05_Order)

#Filter by Reference
#df.primate <- df.primate %>% mutate(found_ref = grepl("2655", References))
#df.primate <- df.primate %>% filter(found_ref == TRUE)

counts <- df.primate %>% group_by(df.primate$MSW05_Genus) %>% count() %>% arrange(desc(n))
#print(counts)

df.primate

#Dags

primate.dag <- dagitty("dag {
  
  Order -> Family
  Family -> Genus
  Genus -> Species
  
  Order -> Forearm_Length
  Order -> Head_Body_Length
  Order -> Body_Mass
  
  Family -> Forearm_Length
  Family -> Head_Body_Length
  Family -> Body_Mass
  
  Genus -> Forearm_Length
  Genus -> Head_Body_Length
  Genus -> Body_Mass
  
  Species -> Forearm_Length
  Species -> Head_Body_Length
  Species -> Body_Mass

  Forearm_Length -> Body_Mass
  Head_Body_Length -> Body_Mass
  
  Forearm_Length <-> Head_Body_Length
       
}")

plot(graphLayout (primate.dag))

adjustmentSets(primate.dag, exposure = c("Order","Family", "Genus", "Species", "Head_Body_Length"), outcome = "Body_Mass", effect="total")
impliedConditionalIndependencies(primate.dag)
## Family _||_ Species | Genus
## Genus _||_ Order | Family
## Order _||_ Species | Genus
## Order _||_ Species | Family

This dag was narrowed down due to the lack of data in each species and no measurments of forearm length.

primate.dag <- dagitty("dag {
  
  Order -> Family
  Family -> Genus

  Order -> Head_Body_Length
  Order -> Body_Mass
  
  Family -> Head_Body_Length
  Family -> Body_Mass
  
  Genus -> Head_Body_Length
  Genus -> Body_Mass
  
  Head_Body_Length -> Body_Mass
  
}")

plot(graphLayout (primate.dag))

adjustmentSets(primate.dag, exposure = c("Order","Family", "Genus", "Species", "Head_Body_Length"), outcome = "Body_Mass", effect="total")
##  {}
impliedConditionalIndependencies(primate.dag)
## Genus _||_ Order | Family
#df.primate <- df.primate %>% mutate(scaled_AdultForearmLen_mm = scale(df.primate$X8.1_AdultForearmLen_mm))
df.primate <- df.primate %>% mutate(scaled_Adult_Mass_g = scale(df.primate$X5.1_AdultBodyMass_g))
df.primate <- df.primate %>% mutate(scaled_Adult_Length_g = scale(df.primate$X13.1_AdultHeadBodyLen_mm))

#Visualize the data

ggplot(df.primate, aes(y=X5.1_AdultBodyMass_g, x=X13.1_AdultHeadBodyLen_mm , color=MSW05_Family)) + geom_point() + coord_cartesian(xlim = c(100,900), ylim = c(-5,50000))

ggplot(df.primate, aes(y=X5.1_AdultBodyMass_g, x=X13.1_AdultHeadBodyLen_mm , color=MSW05_Family)) + geom_point()

ggplot(df.primate, aes(y=X5.1_AdultBodyMass_g, x=X13.1_AdultHeadBodyLen_mm , color=MSW05_Genus)) + geom_point() + coord_cartesian(xlim = c(100,900), ylim = c(-5,50000)) +theme(legend.position="none")

ggplot(df.primate, aes(y=X5.1_AdultBodyMass_g, x=X13.1_AdultHeadBodyLen_mm , color=MSW05_Genus)) + geom_point() + theme(legend.position="none")

This visualization shows that there are clusters for each Primate Family as well as each Genus. These will be used to model a hierachical relationship with the head-body lengths in a multilevel regression.

#Priors and Model

get_prior(scaled_Adult_Mass_g ~ X13.1_AdultHeadBodyLen_mm + (1+X13.1_AdultHeadBodyLen_mm|MSW05_Family/MSW05_Genus) , data=df.primate, family = "gaussian")
my.priors =  prior("normal(0,1)", class= "b" ) + prior("normal(0,1)", class = "Intercept") + prior("normal(0,1)", class="sd")
fit.prior <- brm(scaled_Adult_Mass_g ~ scaled_Adult_Length_g + (1+scaled_Adult_Length_g|MSW05_Family/MSW05_Genus), data=df.primate, family = "gaussian", iter=2000, chains = 4, cores = 4, prior = my.priors, sample_prior = "only" )
## Compiling the C++ model
## Start sampling
plot(fit.prior, ask=F)

pp_check(fit.prior)
## Using 10 posterior samples for ppc type 'dens_overlay' by default.

summary(fit.prior)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: scaled_Adult_Mass_g ~ scaled_Adult_Length_g + (1 + scaled_Adult_Length_g | MSW05_Family/MSW05_Genus) 
##    Data: df.primate (Number of observations: 179) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~MSW05_Family (Number of levels: 15) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.80      0.60     0.03     2.17 1.00
## sd(scaled_Adult_Length_g)                0.80      0.61     0.03     2.24 1.00
## cor(Intercept,scaled_Adult_Length_g)    -0.00      0.56    -0.94     0.93 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            3818     2435
## sd(scaled_Adult_Length_g)                4060     2381
## cor(Intercept,scaled_Adult_Length_g)     8141     2279
## 
## ~MSW05_Family:MSW05_Genus (Number of levels: 67) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.81      0.62     0.03     2.28 1.00
## sd(scaled_Adult_Length_g)                0.81      0.59     0.04     2.23 1.00
## cor(Intercept,scaled_Adult_Length_g)    -0.00      0.58    -0.96     0.95 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            3608     1752
## sd(scaled_Adult_Length_g)                3849     2157
## cor(Intercept,scaled_Adult_Length_g)     9277     1898
## 
## Population-Level Effects: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                -0.03      0.97    -1.94     1.81 1.00     9305
## scaled_Adult_Length_g    -0.01      0.99    -1.97     1.94 1.00    10641
##                       Tail_ESS
## Intercept                 2839
## scaled_Adult_Length_g     2748
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    10.88     12.66     0.33    40.61 1.00     5840     1989
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

This prior model starts out assuming there is no relationship between head-body length and mass.

plot(marginal_effects(fit.prior), ask=F)
## Warning: Method 'marginal_effects' is deprecated. Please use
## 'conditional_effects' instead.

#Vespertilionidae   Eptesicus   nasutus

df.pred <- data.frame(scaled_Adult_Length_g =  seq(-2,2,0.1), MSW05_Family="Cercopithecidae", MSW05_Genus="Cercocebus" )
df.pred
pred <- predict(fit.prior, newdata = df.pred)

df.pred.all <-cbind(df.pred, pred)

ggplot(df.pred.all, aes(x=scaled_Adult_Length_g,y=Estimate)) + geom_point()

Fit Model With Data

fit.bones <- brm(scaled_Adult_Mass_g ~ scaled_Adult_Length_g + (1+scaled_Adult_Length_g|MSW05_Family/MSW05_Genus), data=df.primate, family = "gaussian", iter=2000, chains = 4, cores = 4, prior = my.priors)
## Compiling the C++ model
## recompiling to avoid crashing R session
## Start sampling
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(fit.bones)
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: scaled_Adult_Mass_g ~ scaled_Adult_Length_g + (1 + scaled_Adult_Length_g | MSW05_Family/MSW05_Genus) 
##    Data: df.primate (Number of observations: 179) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~MSW05_Family (Number of levels: 15) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.37      0.10     0.21     0.58 1.00
## sd(scaled_Adult_Length_g)                0.38      0.08     0.25     0.57 1.00
## cor(Intercept,scaled_Adult_Length_g)     0.98      0.04     0.87     1.00 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            1226     1726
## sd(scaled_Adult_Length_g)                1591     1979
## cor(Intercept,scaled_Adult_Length_g)     1253     1339
## 
## ~MSW05_Family:MSW05_Genus (Number of levels: 67) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.10      0.02     0.06     0.14 1.00
## sd(scaled_Adult_Length_g)                0.09      0.02     0.06     0.14 1.00
## cor(Intercept,scaled_Adult_Length_g)     0.91      0.09     0.64     1.00 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            1558     2329
## sd(scaled_Adult_Length_g)                1595     2321
## cor(Intercept,scaled_Adult_Length_g)     1494     1833
## 
## Population-Level Effects: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                -0.17      0.11    -0.39     0.04 1.00      814
## scaled_Adult_Length_g     0.33      0.11     0.11     0.55 1.00      817
##                       Tail_ESS
## Intercept                 1087
## scaled_Adult_Length_g     1259
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.11      0.01     0.10     0.12 1.00     3349     2930
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit.bones)

pp_check(fit.bones)
## Using 10 posterior samples for ppc type 'dens_overlay' by default.

pp_check(fit.bones, type="intervals")
## Using all posterior samples for ppc type 'intervals' by default.

This assumes the total regression, so we explore each regression by grouped clusters.

plot(marginal_effects(fit.bones), ask=F, points=TRUE)
## Warning: Method 'marginal_effects' is deprecated. Please use
## 'conditional_effects' instead.

Here, we see each multilevel regression by family. As we know, there

Model Visualization

coef(fit.bones)$MSW05_Family
## , , Intercept
## 
##                    Estimate  Est.Error       Q2.5       Q97.5
## Aotidae         -0.37695315 0.21457876 -0.7966475  0.05137255
## Atelidae        -0.02754708 0.05499033 -0.1322643  0.08533033
## Cebidae         -0.28592281 0.10098433 -0.4848308 -0.09027649
## Cercopithecidae -0.04878720 0.03132087 -0.1080289  0.01668930
## Cheirogaleidae  -0.34761766 0.16432183 -0.6873702 -0.03250054
## Daubentoniidae  -0.22211031 0.18518926 -0.5861527  0.13844094
## Galagidae       -0.36537515 0.15421333 -0.6869757 -0.07050006
## Hominidae        0.91745711 0.20267376  0.4860554  1.30213604
## Hylobatidae     -0.14171765 0.06170918 -0.2577011 -0.01282146
## Indriidae       -0.19026778 0.07725071 -0.3481673 -0.03949486
## Lemuridae       -0.31928781 0.07109031 -0.4618774 -0.18038532
## Lepilemuridae   -0.27611776 0.25220679 -0.7954933  0.20983495
## Lorisidae       -0.27017606 0.21479960 -0.7126782  0.13654730
## Pitheciidae     -0.31036943 0.07022229 -0.4479488 -0.17578198
## Tarsiidae       -0.34849240 0.18409603 -0.7179206  0.01566124
## 
## , , scaled_Adult_Length_g
## 
##                  Estimate  Est.Error         Q2.5     Q97.5
## Aotidae         0.1362992 0.25446620 -0.359433121 0.6720299
## Atelidae        0.4351337 0.07699635  0.269288626 0.5763010
## Cebidae         0.2190944 0.09907545  0.024674468 0.4146046
## Cercopithecidae 0.4316964 0.04052116  0.350499839 0.5105852
## Cheirogaleidae  0.1533433 0.12626433 -0.099440899 0.3971016
## Daubentoniidae  0.2786565 0.22827554 -0.162031344 0.7450809
## Galagidae       0.1389521 0.13449936 -0.125071518 0.4107979
## Hominidae       1.4745848 0.09183447  1.297929357 1.6686765
## Hylobatidae     0.3242058 0.06836725  0.181681223 0.4547250
## Indriidae       0.3060953 0.08397554  0.139098999 0.4755495
## Lemuridae       0.1752057 0.09221433 -0.009475324 0.3546204
## Lepilemuridae   0.2345731 0.25500899 -0.269735495 0.7344529
## Lorisidae       0.2376647 0.20691881 -0.171223875 0.6409735
## Pitheciidae     0.1842543 0.09556931 -0.009304264 0.3681418
## Tarsiidae       0.1347348 0.13859287 -0.134224373 0.4100331
intercepts <- coef(fit.bones)$MSW05_Family[1] # Specifying the first column only
slopes <- coef(fit.bones)$MSW05_Family[2] # Specifying the second column only

summary(fit.bones)$MSW05_Family
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
## NULL
ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Family))+
  geom_point(aes(color=MSW05_Family))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Family)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Family, fill=MSW05_Family), method="lm", size=0.75,alpha = 0.2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

#+ # average slope with SE
  coord_cartesian(xlim = c(-2,2), ylim = c(-0.75,1))
## <ggproto object: Class CoordCartesian, Coord, gg>
##     aspect: function
##     backtransform_range: function
##     clip: on
##     default: FALSE
##     distance: function
##     expand: TRUE
##     is_free: function
##     is_linear: function
##     labels: function
##     limits: list
##     modify_scales: function
##     range: function
##     render_axis_h: function
##     render_axis_v: function
##     render_bg: function
##     render_fg: function
##     setup_data: function
##     setup_layout: function
##     setup_panel_guides: function
##     setup_panel_params: function
##     setup_params: function
##     train_panel_guides: function
##     transform: function
##     super:  <ggproto object: Class CoordCartesian, Coord, gg>
ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Family))+
  geom_point(aes(color=MSW05_Family))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Family)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Family, fill=MSW05_Family), method="lm", size=0.75,alpha = 0.2) +
 coord_cartesian(xlim = c(100,900), ylim = c(-5,20000))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Family))+
  geom_point(aes(color=MSW05_Family))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Family)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Family, fill=MSW05_Family), method="lm", size=0.75,alpha = 0.1) +
 coord_cartesian(xlim = c(100,500), ylim = c(-5,5000))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Genus))+
  geom_point(aes(color=MSW05_Genus))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Genus)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Genus, fill=MSW05_Genus), method="lm", size=0.75,alpha = 0.2) + theme(legend.position="none")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Genus))+
  geom_point(aes(color=MSW05_Genus))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Genus)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Genus, fill=MSW05_Genus), method="lm", size=0.75,alpha = 0.2) +
 coord_cartesian(xlim = c(100,900), ylim = c(-5,20000)) + theme(legend.position="none")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Genus))+
  geom_point(aes(color=MSW05_Genus))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Genus)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Genus, fill=MSW05_Genus), method="lm", size=0.75,alpha = 0.1) +
 coord_cartesian(xlim = c(100,500), ylim = c(-5,5000)) + theme(legend.position="none")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

ggplot(df.primate, aes(x=X13.1_AdultHeadBodyLen_mm, y=X5.1_AdultBodyMass_g, group=MSW05_Genus))+
  geom_point(aes(color=MSW05_Genus))+
  stat_smooth(method="lm", se=FALSE, size=0.75, aes(color=MSW05_Genus)) + # slopes for different families
  stat_smooth(aes(color=MSW05_Genus, fill=MSW05_Genus), method="lm", size=0.75,alpha = 0.1) +
 coord_cartesian(xlim = c(100,400), ylim = c(-5,1000)) + theme(legend.position="none")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

#Conclusion

In this project, we evaluated the the relationship of mass and head-body length through multilevel regression, considering Family and Genus for Primates.

We can see that more often than not there is a positive relationship between mass and head-body lengths within groups such as their Family and Genus. We also see how the lack of data makes the model very uncertain, relative to the scale at which each group exists. For example, small Primates have a much smaller error, but relative to their size, it expresses just as much of an issue.

What may help, is extra information like age, sex, location and captive/non-captive, other body measurments such as forearm-lengths, ontop of more data for each species (not the mean). There can be even more variation within a species, and this makes it difficult to say that this method is reliable. Even though this is adult data, the generalized approach must reduce the complexity of the information.

Nonetheless, this multilevel regression model exposed both the general relationship between the taxonomic groups, which can at least narrow the physical thresholds of body-mass in certain groups.